Question 1

sf_uscities = readr::read_csv("D:/LifeInUCSB/Study/GEOG176A/week2/assignment5/geog-176A-labs/data/uscities.csv") %>% 
  filter(city ==("Palo")) %>% 
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
  st_transform(5070) %>% 
  st_buffer(5000) %>% 
  st_bbox() %>% 
  st_as_sfc() %>% 
  st_as_sf()

Question 2-3

st = list.files("D:/LifeInUCSB/Study/GEOG176A/week2/assignment5/geog-176A-labs/data", full.names = TRUE, pattern = "TIF")

s = stack(st) %>% 
  setNames(c(paste0("band", 1:6)))


cropper = sf_uscities %>% 
  st_transform(crs(s))

r = crop(s, cropper)

#Step 3

#The dimensions is 7811 rows and 7681 columns, 6 layers. 
#The CRS is WGS84.
#The cell resolution is x=30, y=30 (meters). 

#Step 4
#The dimensions is 340 rows and 346 columns, 6 layers. 
#The CRS is WGS84.
#The cell resolution is x=30, y=30 (meters). 
#R-G-B (natural color)
par(mfrow = c(1,2))
plotRGB(r, r = 4, g = 3, b = 2,stretch = "lin")
plotRGB(r, r = 4, g = 3, b = 2,stretch = "hist")

#NIR-R-G (fa) (color infared)
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 4, b = 3,stretch = "lin")
plotRGB(r, r = 5, g = 4, b = 3,stretch = "hist")

#NIR-SWIR1-R (false color water focus)
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 6, b = 4,stretch = "lin")
plotRGB(r, r = 5, g = 6, b = 4,stretch = "hist")

#My choice
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 7, b = 1,stretch = "lin")
plotRGB(r, r = 5, g = 7, b = 1,stretch = "hist")

#Colored stretch is a method to process different feature in map to keep a visual track. This method adjusts the brightness of a certain wave band to improve its identification and avoid confusion with other types of data.

Question 4

Q4.1

ndvi = (r$band5- r$band4) / (r$band5 + r$band4)

ndwi = (r$band3 - r$band5) / (r$band3 + r$band5)

mndwi = (r$band3 - r$band6) / (r$band3 + r$band6)

wri = (r$band3 + r$band4) / (r$band5 + r$band6)

swi = 1 / (sqrt(r$band2 - r$band6))

stack = stack(ndvi, ndwi, mndwi, wri, swi) %>% 
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
palette = colorRampPalette(c("blue","white","red"))

plot(stack, col = palette(256))

#description
#These pictures show the same distinction between floods near river and land, but with different coloration.
#Flood area in NDWI,MNDWI and WRI plots are depicted in red. In NDVI and SWI plots, floods are drawed in blue.

Q4.2

thresholding1 = function(x){ifelse(x <= 0, 1, NA)}
thresholding2 = function(x){ifelse(x >= 0, 1, NA)}
thresholding3 = function(x){ifelse(x >= 1, 1, NA)}
thresholding4 = function(x){ifelse(x <= 5, 1, NA)}

flood1 = calc(ndvi,thresholding1)
flood2 = calc(ndwi, thresholding2)
flood3 = calc(mndwi, thresholding2)
flood4 = calc(wri, thresholding3)
flood5 = calc(swi, thresholding4)


floodstack = stack(flood1, flood2, flood3, flood4, flood5) %>%
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
plot(floodstack, col = "blue")

Question 5

set.seed(09072020)

value = getValues(r)

dim(value)
## [1] 117640      6
value = na.omit(value)

#Q5.2
#The dataset has 117640 rows and 6 columns.

scvalue = scale(value)


kmeans12 = kmeans(scvalue, centers = 12, iter.max = 100)
kmeans_raster = r$band1
values(kmeans_raster) = kmeans12$cluster
plot(kmeans_raster)

plot(kmeans_raster)

table = table(values(flood1), values(kmeans_raster))

func_kmeans = function(x){ifelse(x == which.max(table), 1, NA)}
cal_kmeans = calc(kmeans_raster, func_kmeans)
floodstack2 = addLayer(floodstack, cal_kmeans)
plot(floodstack2, col = "blue")

Question 6

6.1

kabletable = cellStats(floodstack, sum)
knitr::kable(kabletable, caption = "Total cells number of each plot", col.names = ("Number"))
Total cells number of each plot
Number
NDVI 6666
NDWI 7212
MNDWI 11939
WRI 8469
SWI 15201 
areakable = kabletable * 900
knitr::kable(areakable, caption = "total area of the flooded cells/ m^2", col.names = c("Area"))
total area of the flooded cells/ m^2
Area
NDVI 5999400
NDWI 6490800
MNDWI 10745100
WRI 7622100
SWI 13680900

6.2

floodblues = calc(floodstack2, fun = sum) 
plot(floodblues, col =blues9)

values(floodblues) <- 
  ifelse(values(floodblues)==0,NA, values(floodblues))
mapview(floodblues)
#The pixelation of maps has changed the calculation method of plot cells.
#The reason why some cells are not even is the inaccurate projection after positioning using geographic coordinates. Some cells are segmented differently, resulting in some plot values that are not even.

Extra credit

floodbuilding = st_point(c(-91.78945, 42.06308)) %>% 
  st_sfc(crs = 4326) %>% 
  st_transform(crs(floodstack2)) %>% 
  st_sf()
print("-91.78945, 42.06308")
## [1] "-91.78945, 42.06308"
raster::extract(floodstack2, floodbuilding)
##      NDVI NDWI MNDWI WRI SWI layer
## [1,]    1    1     1   1   1     1

All of the six maps captured the flood in this location.